Quantum Monte Carlo study of quasi-one-dimensional Bose gases 
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We study the behavior of quasi-one-dimensional (quasi- Id) Bose gases by Monte Carlo techniques, 
i.e., by the variational Monte Carlo, the diffusion Monte Carlo, and the fixed-node diffusion Monte 
Carlo technique. Our calculations confirm and extend our results of an earlier study [Astrakharchik 
et al, cond-mat/0308585 . We find that a quasi- fd Bose gas i) is well described by a Id model 
Hamiltonian with contact interactions and renormalized coupling constant; ii) reaches the Tonks- 
Girardeau regime for a critical value of the 3d scattering length a$d', iii) enters a unitary regime for 
|«3d] — * °° 5 where the properties of the gas are independent of azd and are similar to those of a Id 
gas of hard-rods; and iv) becomes unstable against cluster formation for a critical value of the Id 
gas parameter. The accuracy and implications of our results are discussed in detail. 

PACS numbers: 



INTRODUCTION 



Quasi-ld Bose gases have been realized in highly- 
elongated traps by tightly confining the transverse mo- 
tion of the atoms to their zero point oscillations As 
compared to the 3d case, the role of quantum fluctua- 
tions is enhanced in Id and these systems are predicted 
to exhibit peculiar properties, which cannot be described 
using traditional mean-field theories, but require more 
advanced many-body approaches. Particularly intrigu- 
ing is the strong coupling regime, where, due to repul- 
sion between particles, the quasi-ld Bose gas behaves as 
if it consisted of fictitious spinless fermions. This regime, 
the so called Tonks-Girardeau (TG) gas HHHIl, has 
not been achieved yet, but is one of the main focus ar- 
eas of present experimental investigations in this field 6] . 
An interesting possibility to approach the strongly cor- 
related TG regime is provided by magnetic field induced 
atom-atom Feshbach resonances 0>Q. By utilizing this 
technique one can tune the 3d s-wave scattering length 
a^d, and hence the strength of atom-atom interactions, 
to essentially any value, including zero and ±oo. 

Degenerate quantum gases near a Feshbach resonance 
have recently received a great deal of interest both exper- 
imentally and theoretically. At resonance (\a3d\ oo) 
the 3d scattering cross-section a is fixed by the unitary 
condition, a = 4ir/k 2 , where k is the relative wave vec- 
tor of the two atoms. In this regime it is predicted that, 
if the range ro of the atom-atom interaction potential 
is much smaller than the average interparticle distance, 
the behavior of the gas is universal, i.e., independent of 
the details of the interatomic potential and independent 
of the actual value of a^d 03 • This is known as the 
unitary regime [Tl|, In the case of 3d Bose gases, this 
unitary regime can most likely not be realized in experi- 
ments since three-body recombination is expected to set 
in when a^d becomes comparable to the average interpar- 
ticle distance. Three-body recombination leads to cluster 



formation and hence makes the gas-like state unstable. 
The situation is different for Fermi gases, for which the 
unitary regime has already been reached experimentally 
IK ■ In this case, the Fermi pressure stabilizes the system 
even for large |a3<j|. 

In quasi-ld geometries a new length scale becomes rel- 
evant, namely, the oscillator length a p — \Jfij (mw p ) of 
the tightly-confined transverse motion, where m is the 
mass of the atoms and uj p is the angular frequency of the 
harmonic trapping potential. For \a^d\ S> a,p, the gas is 
expected to exhibit a universal behavior if the range ro of 
the atom-atom interaction potential is much smaller than 
a p and the mean interparticle distance is much larger 
than a p. It has been predicted that three-body recom- 
bination processes are suppressed for strongly interact- 
ing Id Bose gases These studies raise the question 
whether the unitary regime can be reached in Bose gases 
confined in highly-elongated traps, that is, whether the 
quasi-ld bosonic gas-like state is stable against cluster 
formation as a^d — > ±00. 

This paper, which is an extension of an earlier study 
[l3| . investigates the properties of a quasi-ld Bose gas at 
zero temperature over a wide range of values of the 3d 
scattering length a^d using quantum Monte Carlo (MC) 
techniques. We find that the system i) is well described 
by a Id model Hamiltonian with contact interactions and 
renormalized coupling constant 4] for any value of a^d', 
ii) reaches the regime of a TG gas for a critical positive 
value of the 3d scattering length a^d] iii) enters a unitary 
regime for large values of |a3d|, that is, for |<i3d| — * 00, 
where the properties of the quasi-ld Bose gas become 
independent of the actual value of a^d and are similar to 
those of a Id gas of hard-rods; and iv) becomes unstable 
against cluster formation for a critical value of the Id gas 
parameter, or equivalently, for a critical negative value of 
the 3d scattering length a^d- 

The structure of this paper is as follows. Section II dis- 
cusses the energetics of two bosons in quasi-ld harmonic 
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traps. We review the mapping of the 3d Hamiltonian to 
a Id model Hamiltonian with contact interactions and 
renormalized coupling constant Q. The eigenenergies 
of the system are calculated by exact diagonalization of 
both the 3d and the Id Hamiltonian. We use these re- 
sults for two particles to benchmark our quantum MC 
calculations presented in Sec. IV. Section III discusses the 
relation between the 3d and the Id Hamiltonian for N 
bosons under quasi-ld confinement. Section IV discusses 
the quantum MC techniques used in the present study: 
variational, diffusion and fixed-node diffusion MC. The 
trial wave functions used for the variational estimates 
and for importance sampling are discussed in detail. Sec- 
tion V presents our MC results for iV = 2 and N = 10 
atoms in highly-elongated harmonic traps over a wide 
range of values of the 3d scattering length 03^. A com- 
parison of the energetics of the lowest-lying gas- like state 
for the 3d and the Id Hamiltonian is carried out. In the 
N = 2 case, we additionally compare with the essentially 
exact results presented in Sec. II. In the N — 10 case, we 
additionally compare with the energy of the lowest-lying 
gas-like state of the Id Hamiltonian calculated using the 
local density approximation (LDA) . Section VI discusses 
the stability of the lowest-lying gas-like state against clus- 
ter formation when a 3 d is negative using the variational 
Monte Carlo (VMC) method. We provide a quantitative 
estimate of the criticality condition. Finally, Sec. VII 
draws our conclusions. 
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which reproduce the low energy scattering solutions of 
the full 3d Hamiltonian, Eq. QJ. Here, £(•) denotes the 
Riemann zeta function, £(1/2) = —1.4604. Alternatively, 
gid can be expressed through the effective Id scattering 
length aid B, 



9id 
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Olshanii's result shows that the waveguide gives rise to an 
effective interaction, parameterized by the coupling con- 
stant gid, which can be tuned to any strength by chang- 
ing the ratio between the 3d scattering length 03^ and 
the transverse oscillator length a p . 

Note that a recent K-matrix treatment results in 
an explicit energy dependence of the effective Id coupling 
constant, 
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II. TWO BOSONS UNDER 
QUASI-ONE-DIMENSIONAL CONFINEMENT 

Consider two interacting mass m bosons with position 
vectors f*i and f%, where r% = (x* in a waveguide 
with harmonic confinement in the radial direction. If wc 
introduce the center of mass coordinate R and the rela- 
tive coordinate r = r% — r±, the problem separates. Since 
the solution to the center of mass Hamiltonian is given 
readily, we only consider the internal Hamiltonian Htff, 
which can be conveniently written in terms of cylindrical 
coordinates r — {p, <f>, z), 
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(1) 



where fi denotes the reduced two-body mass, fi — m/2, 
and V{r) denotes the full 3d atom-atom interaction po- 
tential. 

Considering a regularized zero-range pseudo-potential 



V(f) = 2irh z a 3d /fi6(r)£r, 
scattering length, Olshanii 
Hamiltonian, 
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derives an effective Id 
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where £(•, •) denotes the Hurwitz zeta function and E™ d 
the internal energy of the two boson system. For scat- 
tering between two bosons with minimal internal energy, 
EV^ = hu! p , the Hurwitz zeta function in Eq. ijfjj re- 
duces to the Riemann zeta function, since C(0-5, 1) = 
£(0.5). With the additional assumption that the energy- 
dependence of the 3d scattering length a 3 d can be ne- 
glected, Eq. © reduces to Eq. J2J). For low-energy scat- 
tering the energy-independent effective Id coupling con- 
stant, Eq. ||3J|, is expected to provide a good description 
of the system, and we hence use it throughout this paper. 

The renormalized coupling constant, Eq. (J3J, can be 
compared with the unrenormalized coupling constant 

g°u, 



2?i a 3 d 

gid — 2 ' 



(7) 



which is obtained by averaging the 3d coupling constant 
g 3 d = 4:irh 2 a 3 d/m over the transverse oscillator ground 
state (see, e.g., Ref. p|). Figure Q shows the unrenor- 
malized coupling constant g° d [dashed line, Eq. Q] to- 
gether with the renormalized coupling constant [solid 
line, Eq. |(3J)]. For \a 3 d\ -C a p , the renormalized cou- 
pling constant gid is nearly identical to the unrenor- 
malized coupling constant g® d . For large |o3d|, how- 
ever, the confinement induced renormalization becomes 
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FIG. 1: One-dimensional coupling constants gid [Eq. 
solid line] and g\ d [Eq. {0, dashed line] as a function of the 
3d scattering length a^d/ip- The vertical arrow indicates the 
value of for which gi d diverges, a^ d = 0.9684a p . The hor- 
izontal arrow indicates the asymptotic value of gid as \a$d\ — * 
°°, gid = — 1.9368a p 7kd p . Inset: One-dimensional scattering 
length aid, Eq. (HJ, as a function of a^d/a P - The vertical ar- 
row indicates the value of a%d for which aid goes through zero, 
a M = 0.9684a p . The horizontal arrow indicates the asymp- 
totic value of aid as \a,3d\ — > oo, aid = 1.0326a p . The angular 
frequency u) p determines the frequency v p , lu p — 2-kv p (also, 

TtWn — hVn). 



important, and the effective Id coupling constant gid, 
Eq. J2J), has to be used. At the critical value a^ d = 
0.9684a p (indicated by a vertical arrow in Fig.^l, gi d di- 
verges. For azd — > ±oo, gid reaches an asymptotic value, 
gid = — 1.9368<Xp7kJp (indicated by a horizontal arrow in 
Fig. [TJ. Finally, gid is negative for all negative 3d scat- 
tering lengths. The inset of Fig.^shows the effective Id 
scattering length aid, Eq. iJSJ, as a function of a 3d - For 
small positive 03d, aid is negative and it changes sign 



at a 3d = a c 3d (aid = for a 3d 



i% d )- Moreover, aid 
reaches, just as gid, an asymptotic value for \a 3 d\ — ► 00, 
aid = 1.0326a p (indicated by a horizontal arrow in the 
inset of Fig. [TJ. The renormalized Id scattering length 
aid is positive for negative a 3 d, and approaches +00 as 
0.3d —> 0~. Figure suggests that tuning the 3d scatter- 
ing length a 3 d to large values allows a universal quasi-ld 
regime, where gid and aid are independent of a 3 d, to be 
entered. 

The effective coupling constant gid, Eq. J3J), has been 
derived for a wave guide geometry, that is, with no axial 
confinement. However, it also describes the scattering 
between two bosons confined to other quasi-ld geome- 
tries. Consider, e.g., a Bose gas under harmonic confine- 
ment. If the confinement in the axial direction is weak 
compared to that of the radial direction, the scattering 
properties of each atom pair are expected to be described 
accurately by the effective coupling constant gid and the 



effective scattering length aid- 

The internal motion of two bosons under highly- 
elongated confinement can be described by the following 
3d Hamiltonian 



H 3d = 



1 



(8) 



where u z denotes the angular frequency in the longitu- 
dinal direction, u> z — Xlu p (A denotes the aspect ratio, 
A <C 1). The eigenenergies E^ and eigenfunctions tp^f 
of this Hamiltonian satisfy the Schrodinger equation, 



H^f d t {p,z) = E i ^i n d t {p,z). 



The corresponding Id Hamiltonian reads 



(9) 



h 2 d 2 



Hid = ~Y^ +9ld5{z) + 2^ zZ + ^ (10) 

The Id eigenenergies E™ d of the time-independent 
Schrodinger equation, 



(11) 



can be determined semi-analytically by solving the tran- 
scendental equation |l5j. 

gid = 2V2 r ^ Xz + + J^ tan(7rx z ) huj z a z , (12) 

self consistently for Xz (for a given gid)- In the above 
equation, Xz is an effective (possibly non-integer) quan- 
tum number, which determines the energy E z , 



Xz 



E z 

ihw z 



(13) 



The energy E z , in turn, determines the internal Id 

Id 



eigenenergies E" d , 



Elf = \E Z + hw p . (14) 

In Eq. (|12|l . a z denotes the characteristic oscillator length 
in the axial direction, a z = y/h/(mu} z ). 

To compare the eigenenergies E™ d and E™ d , we use, 
for the 3d atom-atom interaction potential V(r), a short- 
range (SR) model potential V (r) that can support 
two-body bound states, 



V SR (r) 



-V 



cosh (r/fo) 



(15) 



In the above equation, Vo denotes the well depth, and tq 
the range of the potential. In our calculations, tq is fixed 
at a value much smaller than the transverse oscillator 
length, ro = 0.la p . To simulate the behavior of a 3 d near 
a field-dependent Feshbach resonance, we vary the well 
depth Vq, and consequently, the scattering length 03^- 
It has been shown that such a model describes many 
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FIG. 2: Three-dimensional s-wave scattering length a^d as 
a function of the well depth Vb for the short-range model 
potential V s R , Eq. 1151 . Each time the 3d scattering length 
diverges a new two-body s-wave bound state is created. Inset: 
Enlargement of the well depth region used in our calculations. 



atom-atom scattering properties near a Feshbach reso- 
nance properly [16|, Figure El shows the dependence of 
the 3d scattering length 03^ on Vq. Importantly, a^d di- 
verges for particular values of the well depth Vb. At each 
of these divergencies, a new two-body s-wave bound state 
is created. The inset of Fig. [21 shows the range of well 
depths Vq used in our calculations. 

To benchmark our MC calculations (see Sees. I1VI and 
Ivj) . we solve the 3d Schrodinger equation, Eq. 
with A = 0.01 for various well depths Vb using a two- 
dimensional B- spline basis in p and z. Figure |3| shows 
the resulting 3d eigenenergies Etff (diamonds) as a func- 
tion of the 3d scattering length a$d- We distinguish be- 
tween two sets of states: i) States with EJ£f > hcu p are 
referred to as gas-like states; their behavior is, to a good 
approximation, characterized by the 3d scattering length 
03d, and is hence independent of the detailed shape of 
the atom-atom potential. The energies of the gas-like 
states are shown in Fig. |3a) . ii) States with E^ < fiWp 
are referred to as molecular-like bound states; their be- 
havior depends on the detailed shape of the atom-atom 
potential. The energies of these bound states are shown 
in Fig.|3Jb). The well depth Vb of the short-range inter- 
action potential V SR is chosen such that V SR supports 

in the absence of the confining potential — no s-wave 
bound state for a^d < 0, and one s-wave bound state 
for azd > 0. Figure Elb) shows that the bound state re- 
mains bound for 1 0.3^ | — * 00 and for negative a%d if tight 
radial confinement is present. In addition, a dashed line 



shows the 3d binding energy, — % /(ma^ d ), which accu- 
rately describes the highest-lying molecular bound state 
in the absence of any external confinement if 03^ is much 
larger than the range ro of the potential V SR . 

The B-spline basis calculations yield not only the in- 



ternal 3d eigenenergies E^, but also the corresponding 
wave functions ^3™;*. The nodal surface of the lowest- 
lying gas-like state, which is to a good approximation 
an ellipse in the pz-plane, is a crucial ingredient of our 
many-body calculations. Sections IIVI and [V] discuss in 
detail how this nodal surface is used to parametrize our 
trial wave function entering the MC calculations. 

To compare the energy spectrum for iV = 2 of the 
effective Id Hamiltonian with that of the 3d Hamiltonian, 
Fig. additionally shows the Id eigenenergies Ef^ (solid 
lines) obtained by solving the Schrodinger equation for 
-ffJ2*> Eq- P0|l . semi- analytically [using the renormalized 
coupling constant t/id, Eq. ©]. Figure|3fa) demonstrates 
excellent agreement between the 3d and the Id internal 
energies for all states with gas-like character. For positive 
(13d, the effective Id Hamiltonian fails to reproduce the 
energy spectrum of the molecular-like bound states of 
the 3d Hamiltonian accurately [see Fig. 01b), and also 
Ref. [Bill]. 

Our main focus is in the lowest-lying energy level 
with gas-like character. This energy branch is shown 
in Fig. |3{c) on an enlarged scale. A horizontal dashed 
line shows the lowest internal 3d eigenenergy for two 
non-interacting spin-polarized fermions (where the anti- 
symmetry of the wave function enters in the z coordi- 
nate). Our numerical calculations confirm |17| that for 
o-3d = a-zd (ffirf ~~ * 00 ) the- two boson system behaves 
as if it consisted of two non-interacting spin-polarized 
fermions (TG gas). The energy is larger than 

that of two non- interacting fermions for a% d > a% d , and 
approaches the first excited state energy of two non- 
interacting bosons for 03^ — > 0~ [indicated by a dotted 
line in Fig. EJa)]. 

For positive gid, the Id Schrodinger equation, Eq. 
does not support molecular- like bound states. Conse- 
quently, the wave function of the lowest-lying gas-like 
state is positive definite everywhere. For negative gid, 
however, one molecular-like two-body bound state exists. 
If aid a z the bound-state wave function is approxi- 
mately given by the eigenstate ipid °f the Id Hamiltonian 
without confinement, Eq. @, 



$3«(z) = exp --L-L 



E 

aid 



with eigenenergy E^Q*, 



TlUJr, 



in a 



(16) 



(17) 



id 



For the highly-elongated trap with A = 0.01 shown 
in Fig. OUb) and positive a\ d the above binding en- 
ergy nearly coincides with the exact eigenenergy of the 
molecular-like bound state obtained from the solution of 
the transcendental equation (|12[) (solid line). The two- 
body binding energy, Eq. i|17fl . is largest for aid — * + 
(gid — > — 00)] m this case, the molecular-like bound state 
wave function is tightly-localized around z — 0, where 
z = Z2 — z±. Consider a system with aid *C a z . For 
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FIG. 3: Internal eigenenergies E mt as a function of the 
3d scattering length a^d/ap for two bosons under highly- 
elongated confinement with A = 0.01. (a) 3d s-wave eigenen- 
ergies -KS* (diamonds) of gas-like states obtained using the 
short-range model potential V s R , Eq. 1)15^ . in a B-spline basis 
set calculation together with internal Id eigenenergies SJJ 
(solid lines). Excellent agreement between the 3d and Id 
energies is found. Horizontal dotted lines show the lowest 
internal eigenenergies for two non-interacting spin-polarized 
bosons, while horizontal dashed lines show those for two non- 
interacting spin-polarized fermions (indicated respectively by 
"B" and "F" on the right hand side), (b) Binding energy of 
molecular-like bound states. In addition to the 3d and Id 
energies [diamonds and solid lines, respectively; see (a)], a 
dashed line shows the 3d binding energy —fi 2 /{jna\ d ). (c) 
Enlargement of the lowest-lying gas-like state. In addition 
to the 3d and Id energies shown in (a), asterisks show the 
3d energies for the interaction potential V SR calculated using 
the FN-DMC technique, and squares the Id energies for the 
contact interaction potential calculated using the (FN-)DMC 
technique. The statistical uncertainty of the (FN-)DMC ener- 
gies is smaller than the symbol size. Good agreement between 
the (FN-)DMC energies (asterisks and squares) and the en- 
ergies calculated by alternative means (diamonds and solid 
lines) is found. 



negative gid (positive aid), the nodes along the relative 
coordinate z of the lowest- lying gas-like wave function (in 
this case, the first excited state) are then approximately 
given by ±a 1( ;. Thus, imposing the boundary condition 
^id = at |z| = a\d and restricting the configuration 
space to z > aid allows one to obtain an approximation 
to the eigenenergy of the first excited eigen state. Fur- 
thermore, imposing the boundary condition ip^ = at 
z = aid is identical to solving the Id Schrodinger equa- 
tion for a hard-rod interaction potential V (z), 



V HH {z) = 



oo 




for z < aid 
for z > aid 



(18) 



For N — 2, asterisks in Fig. EJc) show the fixed-node 
diffusion Monte Carlo (FN-DMC) results obtained using 
the above fictitious hard-rod potential (see Sec. IV Ap . 
Good agreement is found with the exact Id eigenenergies 
obtained from Eqs. l(T2" )l -lfT3 |l . For N > 2 bosons, our 
Id FN-DMC algorithm and our usage of the hard-rod 
equation of state both take advantage of a reduction of 
configuration space similar to that discussed here for two 
bosons (see Sees. IIIII and II V|) . 



III. N BOSONS UNDER 
QUASI-ONE-DIMENSIONAL CONFINEMENT 

For tightly-confined trapped gases the Id regime is 
reached if the transverse motion of the atoms is frozen, 
with all the particles occupying the ground state of the 
transverse harmonic oscillator. At zero temperature, this 
condition requires that the energy per particle is domi- 
nated by the trapping potential, E/N = Huj p + e, where 
the excess energy e is much smaller than the separation 
between levels in the transverse direction, e <C huj p . In 
the following we consider situations where the Bose gas is 
in the Id regime for any value of the 3d scattering length 
a^d- For a fixed trap anisotropy parameter A and a fixed 
number of particles the above requirement is satisfied 
if A^A < 1. For A = 0.01 and N = IQ (as considered in 
Sec. IVBfl this condition is fulfilled. 

To compare the 3d and Id energetics of a Bose gas, we 
consider the 3d and Id Hamiltonian describing A" spin- 
polarized bosons, 



JY 



1 



—W% + -m{wy i+ ^) 



N 



and 

N 

Hid = J2 

2=1 



z_____ 

2m dzf 



± 2 2 

— muj^Z; 
2 z 1 



l<3 



N 

i<j 

Nhuop (20) 
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respectively. The corresponding eigenenergies and eigen- 
functions are given by solving the Schrodinger equations, 



H 3 di>3d(n, r N ) = E 3d ip 3 d(n,- ■ ■ , r N ) 



and 



(21) 



The corresponding gas density is given by the sum of 
squares of single-particle wave functions, 



1 



N-l 



^4l^ 2 (^M) ex P[-(^) 2 ] (25) 



H u ^id{zi, ■ ■ ■ ,z N ) = E ld tp ld (z 1 , ■■■ , z N ), (22) 

respectively. In contrast to Sec. |H] here we do not 
separate out the center of mass motion since the MC 
calculations used to solve the 3d and Id many-body 
Schrodinger equations can be most conveniently imple- 
mented in Cartesian coordinate space (see SecEJ. In 
the following, we refer to eigenstates of the confined Bose 
gas with energy greater than Nhtu p as gas-like states, and 
to those with energy smaller than Nhcu p as cluster-like 
bound states. 

Section lY Bl compares the energetics of the lowest-lying 
gas- like state of the 3d Schrodinger equation, Eq. (|21|l . 
obtained using the short-range potential V s R , Eq. (|15fl . 
with that obtained using the hard-sphere potential V HS , 



V Hts (r) 







for r < a 3d 
for r > a 3d 



(23) 



For V HS , the s-wave scattering length a 3d coincides with 
the range of the potential. For V SR , in contrast, Tq de- 
termines the range of the potential, while the scattering 
length a 3d is determined by r$ and Vq. For a 3d <C a p , 
both potentials give nearly identical results for the en- 
ergetics of the lowest-lying gas-like state, which depend 
to a good approximation only on the value of a 3d . For 
a 3 d ^ Op, instead, deviations due to the different effective 
ranges become visible and only V yields results, which 
do not depend on the short-range details of the potential 
and which are compatible with a Id contact potential. 

Section IV Bl also discusses the energetics of the Id 
Hamiltonian, Eq. (|20|l . For small |<?id|, the energetics 
of the many-body Id Hamiltonian are described well by 
a Id mean-field equation with non-linearity. For nega- 
tive gi d , the mean-field framework describes, for exam- 
ple, bright solitons |20l l2l| . which have been observed 
experimentally J^]. For large |<7id|, in contrast, the sys- 
tem is highly-correlated, and any mean-field treatment 
will fail. Instead, a many-body description that incorpo- 
rates higher order correlations has to be used. In par- 
ticular, the limit \gi d \ — > 00 corresponds to the strongly- 
interacting TG regime. 

For infinitely strong particle interactions {\gi d \ — > 00), 
Girardeau shows [3j, using the equivalence between the 
Id 5-function potential and a "Id hard-point potential" , 
that the energy spectrum of the Id Bose gas coincides 
with that of N non- interacting spin-polarized fermions. 
The lowest eigenenergy per particle of the Id Bose gas, 
Eq. JUJ), is, in the TG limit, given by 



E TG 
N 



1 ] TllOr, 



(24) 



with the normalization f_ nJ G (z)dz = N. In Eq. i|25[l . 
the Hk denote Hcrmitc polynomials, and z denotes the 
distance measured from the center of the trap. For large 
numbers of atoms, the density expression in Eq. (|25|l can 
be calculated using the LDA [23| . 



^^d G (z) 



1 - 



2Nal 



1/2 



(26) 



The above result cannot reproduce the oscillatory behav- 
ior of the exact density, Eq. Q25p. but it does describe the 
overall behavior properly (see Sec. IVI|) . 

To characterize the inhomogeneous Id Bose gas fur- 
ther, we consider the many-body Hamiltonian of the ho- 
mogeneous Id Bose gas, 



H ld = 



N 

^ 2m dzf 

i—l fc 



N 



+ .9idV S(z i:j ) + NJkJ p . (27) 



i<j 



[By introducing the energy shift Nhuo P: our classification 
of gas-like states and cluster-like bound states introduced 
after Eq. (|2^|l remains valid.] For positive gid, H^ m 
corresponds to the Lieb-Liniger (LL) Hamiltonian. The 
gas-like states of the LL Hamiltonian, including its ther- 
modynamic properties, have been studied in detail |24| . 
The energy per particle of the lowest-lying gas-like state, 
the ground state of the system, is given by 



Etd L M 



N 



2m 



e(7)«id! 



(28) 



where ni d denotes the density of the homogeneous sys- 
tem, and e(7) a function of the dimensionless parameter 
7 = 2/(n ld \aid\). 

We use the known properties of the LL Hamiltonian 
to determine properties of the corresponding inhomoge- 
neous system, Eq. (|20fl . within the LDA. This approxi- 
mation provides a correct description of the trapped gas 
if the size of the atomic cloud is much larger than the 
characteristic length scale a z of the confinement in the 
longitudinal direction (23j. Specifically, consider the lo- 
cal equilibrium condition, 

fi(N) = hujp + fiiocal[nid(z)] + ^mLu^z 2 , (29) 

where /i; OC a/("-id) denotes the chemical potential of the 
homogeneous system with density nid, 



{nid) 



d [n ld E^ d L {n ld )/N] 



dn 



id 



(30) 
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The chemical potential (J,(N), Eq. (|29(l . can be calculated 
using Eq. (|3U|) together with the normalization of the den- 
s ity> J^°oo nid{z)dz = N. Integrating the chemical poten- 
tial n(N) then determines the energy of the lowest-lying 
gas-like state of the inhomogeneous A-particle system 
within the LDA. The LDA treatment is computationally 
less demanding than solving the many-body Schrodinger 
equation, Eq. I|22[l . using MC techniques. By compar- 
ing with our full Id many-body results we establish the 
accuracy of the LDA (see Sec. EH- 

For negative gid, the Hamiltonian given in Eq. i|27|) 
supports cluster-like bound states ("tightly-bound 
droplets"). The grou nd state energy and eigenfunction 
of the system are [25| 



E ld 

N 



6ma 2 d 



{N 2 - 1) + Tilu p 



(31) 



and 



N 



^r(^---^A0 = n cx p 



1<J 



aid 



(32) 



respectively. The eigenstate given by Eq. (|32|l depends 
only on the relative coordinates zy, that is, it is inde- 
pendent of the center of mass of the system. Adding a 
confinement potential [see Eq. (|2(J|l ] with uj z such that 
a z ^> aid leaves the eigenenergy E^ m of this cluster-like 
bound state to a good approximation unchanged, while 
the corresponding wave function becomes localized at the 
center of the trap. This state describes a bright soliton, 
whose energy can also be determined within a mean- 
field framework [HI . An excited state of the many-body 
Id Hamiltonian with confinement corresponds, e.g., to a 
state, where N — 1 particles form a cluster-like bound 
state, i.e., a soliton with N — 1 particles, and where 
one particle approximately occupies the lowest harmonic 
oscillator state, i.e., has gas- like character. Similarly, 
molecular-like bound states can form with fewer parti- 
cles. 

The above discussion implies that the lowest-lying 
gas-like state of the Id Hamiltonian with confinement, 
Eq. (|2*0"|) with negative gid, corresponds to a highly- 
excited state. For dilute Id systems with negative gid, 
the nodal surface of this excited state can be well ap- 
proximated by the following nodal surface: ipid = for 
Zij = aid, where i.j ~ 1, • • • , N and i < j. As in the two- 
body case, the many-body energy can then be calculated 
approximately by restricting the configuration space to 
regions where the wave function is positive. This corre- 
sponds to treating a gas of hard-rods of size aid- In the 
low density limit, we expect that the lowest-lying gas- like 
state of the Id many-body Hamiltonian with gid < is 
well described by a system of hard-rods of size aid- 
In addition to treating the full Id many-body Hamilto- 
nian, we treat the inhomogeneous system with negative 
gid within the LDA. The equation of state of the uniform 



hard-rod gas with density nid is given by Q 



N 



n 2 h 2 n 2 



id 



6m (1 - ni d aid) 2 



TlLOn 



(33) 



We use this energy in the LDA treatment [see Eqs. 
through (J30J) with E{ d L replaced by E^ d R }. The hard- 
rod equation of state treated within the LDA provides a 
good description when gid is negative, but \gid\ not too 
small (see Sees. EB and ED- To gain more insight, we 
determine the expansion for inhomogeneous systems with 
NX <C 1 using the equation of state for the homogeneous 
hard-rod gas, 



E u % % NX 



128V2 
45tt 2 



-aid 



— + ■■■} (34) 



The first term corresponds to the energy per particle in 
the TG regime [see Eq. l(2*I)l]: the other terms can be 
considered as small corrections to the TG energy. In the 
unitary limit, that is, for aid/a p = 1.0326, expression 
(|34|l becomes independent of 03^, and depends only on 
NX. Similarly, the linear density in the center of the 
cloud, z — 0, is to lowest order given by the TG result, 
< d G (0) = V2NX/(TTa p ) [see Eq. ®]. Section Ivll shows 
that the TG density provides a good description of in- 
homogeneous Id Bose gases over a fairly large range of 
negative g M . 



IV. QUANTUM MONTE CARLO METHODS 

This section describes the variational, diffusion and 
fixed- node diffusion MC methods (see, e.g., ^26]) used 
in the present study. These quantum MC techniques 
solve the many-body Schrodinger equation for the ground 
state and for excited states at zero temperature. Simi- 
lar to other MC approaches, these techniques are based 
on stochastic numerical algorithms, which are powerful 
when one is treating systems with many degrees of free- 
dom. 



A. Variational Monte Carlo method 

The VMC method was first introduced in the seminal 
work by McMillan [23] to study the ground state prop- 
erties of liquid 4 He. The VMC technique is based on the 
variational principle, 



(ip T \H\ifr) 
(tp T \ip T ) 



> En 



(35) 



where "0T denotes a variational or trial wave function, 
which is parameterized in terms of a set of variational 
parameters. In Eq. (|35|) . H denotes the Hamiltonian of 
a bosonic N particle system with ground state energy 
Eq and stationary ground state wave function ^q- The 
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evaluation of the high-dimensional integral, Eq. (|35|l . can 
be performed by MC techniques, resulting in the VMC 
energy expectation value. This variational energy expec- 
tation value is an upper bound to the exact ground state 
energy Eq. Importantly, the variational principle also ap- 
plies to excited states. For a trial wave function ipT with 
a given symmetry, the variational estimate provides an 
upper bound to the energy of the lowest excited state of 
the Hamiltonian H with that symmetry. 

Choosing a good functional form for the trial wave 
function tpT is a crucial step of the VMC method (and 
also of the P MC me thod and the FN-DMC method, see 
Sees. IIV Bl and llV C|l . A general ansatz for ipx, which has 
been used successfully to describe a system of N spinless 
bosons under external confinement in either 3d or Id m, 
is a Bijl-Jastrow decomposition into one- and two-body 
correlation factors. The one-body term accounts for the 
confining potential, while the two-body term accounts for 
interactions between particles. 

To describe the lowest-lying gas-like state of H^d 
[Eq. (|19|l ]. we use the following trial wave function 



N 



N 



^ T d (r u ... 1 r N ) = Y[, 



-z?/(2al) p -p*/(2 a l) 



i=l 



(36) 

Here, a z and a p determine the Gaussian width of tp^ 1 
in the longitudinal and transverse direction, respectively. 
These variational parameters are optimized in the course 
of the VMC calculation by minimizing the energy ex- 
pectation value. The two-body correlation factor fair) is 
chosen to closely reproduce the scattering behavior of two 
bosons at low energies. For the hard-sphere potential, we 
use 



hif) = 



i 



o 

a3d/\r\ 



for 
for 



M < a 3d 
\r\ > a zd . 



(37) 



The constraint fa = for r < a^d accounts for the bound- 
ary condition imposed by the hard-sphere potential; it is 
exact even for the many-body system. For the short- 
range potential, we use instead 



h(r) 







1-V 



777 



for 
for 



+ 



z 2 



< 1 
> 1 



(38) 



where a and b denote the lengths of the semi-axes of 
an ellipse. For two particles under highly-elongated con- 
finement, the nodal surface is to a good approximation 
cllipticly shaped (see Sec. IV A)l . Thus, the parameters 
a and b are determined by fitting the elliptical surface 
to the nodal surface obtained by solving the Schrodinger 
equation for N = 2, Eqs. (jHJ) and 0, by performing a 
B-spline basis set calculation. In contrast to V H , the 
constraint f% = in Eq. (|38|l parameterizes the many- 
body nodal surface for V" only approximately. We 
expect that our parameterization leads to an accurate 
description of quasi- Id Bose gases if the average distance 
between particles is much larger than the semi-axes of 



the ellipse. The trial wave functions discussed here in 
the context of our VMC calculations also enter our (FN- 
)DMC calculations (see Sec. lTVBl and llVT^ . 

In the case of the Id Hamiltonian, Eq. (|2UII . we use a 
trial wave function of the form 



JV 



N 



^ T d (z u ...,z N ) = l[e-**/^l[ f 2 (z i - Zj ) , (39) 

where the Gaussian width a z is treated as a variational 
parameter. The two-body correlation factor f 2 {z) is cho- 
sen as 

Mg) = {Mk.(\z\-Z)] for m 

The cut-off length Z is fixed at Z = 500aid, while the 
wave vector k z is chosen such that the boundary condi- 
tion at z = imposed by the (^-function potential is sat- 
isfied: —k z t&n(k z Z) = l/a±d- For negative aid {gid > 0) 
the correlation function, Eq. I|40|l . is positive everywhere. 
For positive aid {gid < 0), in contrast, fviz) changes sign 
at \z\ = aid- The parameterization given by Eq. (|40|l is 
used in our stability analysis performed within a VMC 
framework (see Sec. IVIf) and in our DMC calculations 
for g M > (see Sec. HVB) . To perf orm t he FN-DMC 
calculations for negative gid (see Sec. IIV(J|) . we need to 
construct a trial wave function that is positive definite 
everywhere. In the FN-DMC calculations, we thus use 
an alternative parameterization, which imposes the con- 
straint /2 = for aid < z, 

( for z < aid 

f 2 (z)= { cos[k z {\z\ - Z)] for a ld <z<Z (41) 
1 for z > Z . 



B. Diffusion Monte Carlo method 

The diffusion Monte Carlo (DMC) algorithm solves 
the time-independent Schrodinger equation of a N par- 
ticle system by introducing the imaginary time r = it/h, 
and propagating the function f(lZ,r) — ^t(T^)^(T^,t) 
in imaginary time r. Here, ^(TZ, r) denotes the wave 
function of the system, which we are seeking. The trial 
wave function i/jt{7Z-)-, which can be optimized in a se- 
ries of VMC calculations, enters the DMC calculations 
as input. In 3d, 1Z collectively denotes the position vec- 
tors, 1Z = (r*i, rjv); in Id, it collectively denotes the 
positions, 1Z = (zi, • • • , zjv). 

Using the function f(lZ,r), the time-dependent 
Schrodinger equation in imaginary time can be rewrit- 
ten as 



df(K,' 



= - DV 2 n f(n,T) + DVn[F(n)f(n,T)} 



[E L (TZ) - E ref ]f(n,T) , 



(42) 



where El{TZ) = i/jt(TZ) 1 H4)t(JZ) denotes the local en- 
ergy, F(1Z) = 2iIj t ('JZ)- 1 X7k4>t(TZ) the quantum drift 
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force, and D the diffusion constant, D = h 2 /2m. The 
subscript 7Z of the operator V indicates that the deriva- 
tive has to be taken for every component of 1Z. The 
constant energy shift E re f is introduced to stabilize the 
numerics. The solution of Eq. (|42|l can be written for- 
mally as 

f(Tl', r + At) = J dTZ G(TZ\TZ, Ar)/(ft, r) , (43) 
where the time-dependent Green's function G is given by 

G{n',n,Ar) = (n'\e- AAT \n) . (u) 

Here, A denotes the time evolution operator of Eq. (|42|l . 
— df(TZ, t)/8t = Af(TZ,r). If the short-time approxi- 
mation to the Green's function G(1Z', 1Z, At) is known 
(for sufficiently small At), the asymptotic solution for 
large times t, f(lZ,T — > oo), can be obtained by propa- 
gating / for a large number of time steps At (for more 
details see, e.g., |2ij). For a system of bosons, one can 
show that f(ll,T oo) = ip T (TZ)^ (TZ), where <f is 
the exact ground-state eigenfunction. We calculate the 
eigenenergy Eg by first propagating to large times t, and 
then averaging the local energy El over the distribution 

= J dKMK)Hy (K) = J dlZf(lZ, t — > oo)E L (K) 
j dn^ T {TV)^ a {n) J dTZf(TZ,T ^ oo) 

(45) 

Apart from statistical errors, the DMC method deter- 
mines the energy of the nodeless ground state of a system 
of N bosons essentially exactly. Importantly, the energy 
expectation value calculated by the DMC technique with 
importance sampling, which we use in our study, is in- 
dependent of the detailed shape of i/t as long as ifo is 
positive definite everywhere. 

We use the outlined algorithm to calculate the ground- 
state energy E%d of the 3d Hamiltonian, Eq. (|19|l . for the 
hard-sphere interaction potential, V HS , and that of the 
Id Hamiltonian, Eq. i|2U|) . for positive coupling constant, 
gid > 0. The trial wave function ipx used in these cases 
is given by Eqs. |j2§ and l(37|l. and by Eqs. and (pE0|> - 
respectively. 

The DMC algorithm cannot be used directly to calcu- 
late excited states. If ^ is an excited state, that is, if ^ is 
orthogonal to *S? , the function f(TZ, t) = iPt(TZ)^(T^, t) 
is not positive everywhere in configuration space and 
can hence not be interpreted as a probability density. 
This leads to the so-called sign problem. To neverthe- 
less calculate excited state energies such as the energy 
of the lowest-lying gas-like state of the 3d Hamiltonian, 
Eq. 119|) . with the short-range interatomic potential V s R , 
Eq. JT51), or of the Id Hamiltonian, Eq. J5D|) with g ld < 0, 
we apply the FN-DMC method. 

C. Fixed- node diffusion Monte Carlo method 

The FN-DMC method jlj modifies the DMC method 
to allow approximate treatment of excited states of 



many-body systems. The idea of the FN-DMC method 
is to treat excited states by "enforcing" the positive 
definitcncss of the probability distribution f(lZ 7 r) — 
iPt(TZ)^{TZ,t). The function / is positive definite ev- 
erywhere in configuration space, and can hence be in- 
terpreted as a probability distribution, if -0T and SS* 
change sign together, and thus share the same (high- 
dimensional) nodal surface. To ensure positive definite- 
ness of /, the trial wave function -0T imposes a nodal 
constraint, which is fixed during the calculation. Within 
this constraint, the function / is propagated (following a 
scheme very similar to that outlined in Sec. IIVB|l . and 
reaches an asymptotic distribution for large t, f(lZ,r — ► 
oo) = ip T (TZ)^(TZ). In the FN-DMC method, * is 
an approximation to the exact excited eigenfunction of 
the many-body Schrodinger equation (and not the exact 
eigenfunction as in the DMC method). It can be proven 
that, due to the nodal constraint, the fixed-node energy 
is a variational upper bound to the exact eigenenergy for 
a given symmetry 30] . In particular, if the nodal sur- 
face of ipT were exact, then \I> would be exact. Thus, the 
FN-DMC energy depends crucially on a good parameter- 
ization of the many-body nodal surface. 

Section El reports the energy of the lowest-lying gas- 
like state of the 3d Hamiltonian, Eq. I|19|) . for the short- 
range potential V SR , and that of the Id Hamiltonian, 
Eq. for g xd < calcula ted u sing the FN-DMC 

method. As discussed in Sec. IIV Al the nodal surface 
of the many-body trial wave function ipT is constructed 
from the essentially exact nodal surface of the two-body 
wave function describing the lowest-lying gas-like state. 
In the 3d case, the importance sampling trial wave func- 
tion is given by Eqs. (|3*tj)l and (|3*5|l . In the Id case, the 
importance sampling trial wave function is instead given 
by Eqs. (3HJ! and gTJ. We expect that the FN-DMC ap- 
proach implemented here results in accurate many-body 
energies of dilute quasi- Id Bose gases. 



V. ENERGETICS OF 
QUASI-ONE-DIMENSIONAL BOSE GASES 

Table H] summarizes the techniques used to solve the Id 
and 3d Hamiltonian, respectively. This table is meant to 
guide the reader through our result sections. Sect ion fV Al 
discusses our MC energies for two-particle systems, while 
Sec. IV Bl discusses the energetics for larger systems, cal- 
culated within various frameworks. Finally, Sec. I VII dis- 
cusses the stability of quasi- Id Bose gases. 



A. Two-body system 

Section ITU discusses the calculation of the energy spec- 
trum related to the internal motion of two bosons under 
highly-elongated confinement, Eq. (jHJ, using a B-spline 
basis, and the eigenspectrum related to the internal mo- 
tion of two bosons under Id confinement, Eq. <|1(J|I . using 
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Hamiltonian 


interaction 


technique 


Section 


H:id 


yHS 


DMC 


VB 


H^d 


ySR 


FN-DMC 


VA, VB 


Hid 


gid > 


DMC 


VA, VB 


Hu 


gid > 


LDA, LL 


VB 


Hi d 


gid < 


FN-DMC 


VA, VB 


H ld 


gid > 


LDA, hard-rod 


VB 


Hid 


gid < 


VMC 


VI 



TABLE I: Guide that summarizes the techniques used to 
solve the 3d and Id Hamiltonian, Eqs. 1191 and (121)1 . respec- 
tively. Column 2 specifies the atom-atom interactions of the 
many-body Hamiltonian, column 3 lists the techniques used 
to solve the corresponding many-body Schrodinger equation, 
and column 4 lists the sections that discuss the results ob- 
tained using this approach. 




Eqs. (fl^|l through ifbljl . We now use these essentially ex- 
act eigenenergies to benchmark our (FN-)DMC calcula- 
tions. Toward this end, we solve the 3d Schrodinger equa- 
tion, Eq. (|19fl . and the Id Schrodinger equation, Eq. <|20|) . 
for various interaction strengths for N = 2 and A = 0.01 
using (FN-)DMC techniques. The resulting MC ener- 
gies E^d and E\d include the center of mass energy of 
(1 + \/2)htu p . To compare with the internal eigenener- 
gies discussed in Sec. [H] we subtract the center of mass 
energy from the (FN-)DMC energies. 

For N = 2, the lowest-lying gas- like state of the 3d 
Hamiltonian for the short-range potential V SR is the 
first excited eigenstate. Consequently, we solve the 3d 
Schrodinger equation by the FN-DMC technique using 
the trial wave function ipx given by Eqs. (|36(l and (|38ll . 
Figure 0] shows the elliptical nodal surface of the trial 
wave function ipx (solid lines) together with the essen- 
tially exact nodal surface calculated using a B-spline ba- 
sis set (symbols; see also Sec. |HJ for A = 0.01 and three 
different scattering lengths, a^d/a. p = 1,6, and —4. No- 
tably, the semi-axes a along the p coordinate is larger 
than that along the z coordinate (a/b > 1), "opposing" 
the shape of the confining potential, for which the charac- 
teristic length along the p coordinate is smaller than that 
along the z coordinate (a p /a z < 1). Figure 0] indicates 
good agreement between the essentially exact nodal sur- 
face and the parameterization of the nodal surface by an 
ellipse for a^d/o-p = 1 and 6; some discrepancies become 
apparent for negative a^d- Since the FN-DMC method 
results in the exact energy if the nodal surface of ipT 
coincides with the nodal surface of the exact eigenfunc- 
tion, comparing the FN-DMC energies for two particles 
with those obtained from a B-spline basis set calculation 
provides a direct measure of the quality of the nodal sur- 
face of ipT- Figure Etc) compares the internal 3d energy 
of the lowest-lying gas-like state calculated using a B- 
spline basis (diamonds, see Seem with that calculated 
using the FN-DMC technique (asterisks). The agreement 
between these two sets of energies is — within the sta- 
tistical uncertainty — excellent for all scattering lengths 



FIG. 4: Nodal surface of the trial wave function tpx [solid 
lines, Eq. I|38fl ] together with the essentially exact nodal 
surface calculated using a B-spline basis set (symbols) for 
A = 0.01, N — 2, and three different scattering lengths, 
d3d/a P — 1 (pluses), asd/dp — 6 (asterisks), and a 3d /a p — —4 
(diamonds). The nodal surface is shown as a function of the 
internal coordinates z and p. Good agreement between the 
elliptical nodal surface (solid lines) and the essentially exact 
nodal surface (symbols) is visible for a^d/a-p = 1 and 6. Small 
deviations are visible for asd/cip = —4. 



a,3d considered. We conclude that our parameterization 
of the two-body nodal surface, Eq. I|38|) . is accurate over 
the whole range of interaction strengths 03^ considered. 

We expect that our parameterization of the two-body 
nodal surface is to a good approximation independent 
of the confining potential in z (for small enough A). In 
fact, we expect our nodal surface to closely resemble that 
of the scattering wave function at low scattering energy 
of the 3d wave guide Hamiltonian given by Eq. Q . To 
quantify this statement, Fig.[S]shows the semi-axes a and 
b (pluses and asterisks, respectively) obtained by fitting 
an ellipse, see Eq. I|38l) . to the nodal surface obtained by 
solving the Schrodinger equation for the two-body Hamil- 
tonian, Eq. ©, using a B-splinc basis for various aspect 
ratios A (A = 0.001, •••,!), and fixed scattering length, 
asd = 2a p (similar results are found for other scattering 
lengths). Indeed, the nodal surface for a given a^d/o-p is 
nearly independent of the aspect ratio A for A < 0.01. 
These findings for two particles imply that the param- 
eterization of the nodal surface of ipT used in the FN- 
DMC many-body calculations should be good as long as 
the density along z is small. 

Next, consider the Id Hamiltonian, Eq. I|20|) . for N = 
2. For positive gid, the lowest-lying gas-like state is the 
ground state of the two-body system and we hence use 
the DMC technique [with ipr given by Eqs. i|3D|) and 
(|4*U|l ]; for gid < 0, the lowest-lying gas-like state is the 
first excited state, and we instead use the FN-DMC tech- 
nique [with ipT given by Eqs. (|39|l and i|41|l ]. Figure |3| 
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FIG. 5: Semi-axes a (pluses) and b (asterisks) obtained by 
fitting an ellipse [see Eq. I|38fl ] to the essentially exact nodal 
surface for two bosons under cylindrical confinement, calcu- 
lated using a B-spline basis set as a function of the anisotropy 
parameter A, for a^dja-p = 2. Dotted lines are shown to guide 
the eye. For A < 0.01, the nodal surface is nearly independent 
of the anisotropy parameter A. 



shows the Id energies of the lowest-lying gas-like state 
calculated using Eqs. l(T2|l through ljT4")l (solid line), to- 
gether with those calculated by the (FN-)DMC technique 
(squares). We find excellent agreement between these 
two sets of Id energies. 

The comparison for two bosons between the (FN-)- 
DMC energies and the energies calculated by alternative 
means serves as a stringent test of our MC codes, since 
these codes are implemented such that the number of 
particles enters simply as a parameter. 
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scattering length a3d [a p ] 

FIG. 6: Three-dimensional (FN-)DMC energy per particle, 
Ezd/N - ?uj p , calculated using V HS (diamonds) and V SR 
(asterisks), respectively, together with Id (FN-)DMC energy 
per particle, E 1 d/N — hw p , calculated using gid [squares, 
Eq. J3J] and g® d [pluses, Eq. J7J], respectively, as a function 
°f a 3d [(a) linear scale; (b) logarithmic scale] for N = 10 and 
A = 0.01. The statistical uncertainty of the (FN-)DMC en- 
ergies is smaller than the symbol size. Dotted and solid lines 
show the Id energy per particle calculated within the LDA 
for g® d , Eq. [using the LL equation of state] and for gid, 
Eq. [using the LL equation of state for > 0, and the 
hard-rod equation of state for gid < 0], respectively. A dot- 
ted horizontal line indicates the energy per particle of a non- 
interacting Bose gas, and a dashed horizontal line indicates 
the TG energy per particle. A vertical arrow the position 
where gid, Eq. @, diverges. 



B. N-body system 

This section presents our many-body study, which in- 
vestigates the properties of quasi-ld Bose gases over a 
wide range of scattering lengths a$d- We focus specifi- 
cally on three distinct regimes: i) < a^d < a^ d (gid 
is positive); ii) \asd\ — > oo (gid and aid are independent 
of a^d] unitary regime); and iii) 03^ — * 0~ (aid is large 
and positive; onset of instability). We discuss the ener- 
getics of quasi-ld Bose gases for N = 10. Our results 
presented here support our earlier conclusions, which are 
based on a study conducted for a smaller system, i.e., for 
N = 5 H3. 

For small A (here, A = 0.01), the radial angular fre- 
quency lu p dominates the eigenenergies of the 3d and of 
the Id Schrodinger equation. The shift of the eigenen- 
ergy of the lowest-lying gas-like state as a function of the 
interaction strength is, however, set by the axial angu- 
lar frequency lu z . To emphasize the dependence of the 
eigenenergies of the lowest- lying gas- like state on u z , we 



report the energy per particle subtracting the constant 
offset hujp, that is, we report the quantity E/N — fko p . 

Consider the lowest-lying gas-like state of the 3d 
Schrodinger equation. Figure shows the 3d energy 
per particle, E^d/N — Tluj p , as a function of a$d for 
N = 10 under quasi-ld confinement, A = 0.01, for the 
hard-sphere two-body potential V HS (diamonds) and the 
short-range potential V s R (asterisks). The energies for 
V HS are calculated using the DMC method [with i/jt 
given by Eqs. (|36|l and l|37|) ]. while those for V SR are 
calculated using the FN-DMC method [with ipx given 
by Eqs. I|36() and J3SJ]. For small a^d/a-p, the energies 
for these two two-body potentials agree within the sta- 
tistical uncertainty. For a^d o-p, however, clear dis- 
crepancies are visible. The DMC energies for V SR cross 
the TG energy per particle (indicated by a dashed hor- 
izontal line), E/N — huj p — huj p XN/2, very close to the 
value a^ d = 0.968Aa p [indicated by a vertical arrow in 
Fig. EJb)]j while the energies for V HS cross the TG en- 
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ergy per particle at a somewhat smaller value of a^d- 
For V HS , "fermionization" of a quasi- Id gas has previ- 
ously been investigated in detail [13 ■ The present paper 
goes beyond these previous studies in that we consider a 
short-range potential, whose applicability extends to the 
regime a 3 d > a^ d . 

For 0,3d > a^ d , the energy for the short-range potential 
V SR of the lowest-lying gas-like state increases slowly 
with increasing a%d, and becomes approximately constant 
for large values of |a3d|. The limit \a 3 d\ — * 00 corresponds 
to the unitary regime (see below). Notably, the 3d energy 
behaves smoothly as a 3 d diverges. The 3d energy slowly 
increases further for increasing negative 03d, and changes 
more rapidly as a 3 d — ► 0~. The |a3d| — > 00 regime and 
the a^d —> Q~ regime are discussed in more detail below. 

To compare our results obtained for the 3d Hamilto- 
nian, H%d, with those for the Id Hamiltonian, Hid, we 
also solve the Schrodinger equation for Hid, Eq. (J2DJ, 
for the lowest-lying gas-like state. For positive coupling- 
constants, gid > 0, the lowest-lying gas-like state is the 
many-body ground state, and we hence use the DMC 
method [with ipx given by Eqs. (|3"5|l and 1(10)1]. For 
gid < 0, however, the Id Hamiltonian supports cluster- 
like bound states. In this case, the lowest-lying gas-like 
state corresponds to an excited many-body state, and we 
hence solve the Id Schrodinger equation by the FN-DMC 
method [with ip T given by Eqs. (JSHJ) and {2J]. 

Figure shows the resulting Id energies per particle, 
Eid/N — hup, for the renormalizcd coupling constant gid 
[squares, Eq. ©], and the unrenormalizcd coupling con- 
stant gi d [pluses, Eq. JJJ], respectively. The Id ener- 
gies calculated using the two different coupling constants 
agree well for small a 3d , while clear discrepancies become 
apparent for a%d ^ a^ d . In fact, the Id energies cal- 
culated using the unrenormalized coupling constant g® d 
approach the TG energy (dashed horizontal line) asymp- 
totically for a 3d — * 00, but do not become larger than the 
TG energy. The Id energies calculated using the renor- 
malized Id coupling constant gid agree well with the 3d 
energies calculated using the short-range potential V SR 
(asterisks) up to very large values of the 3d scattering 
length a%d- ln contrast, the Id energies deviate clearly 
from the 3d energies calculated using the hard-sphere po- 
tential V HS (diamonds) at large 03^. 

The Id energies calculated using the renormalized cou- 
pling constant agree with the 3d energies calculated using 
the short-range potential V SR also for [03^ | — > 00, that 
is, in the unitary regime, and for negative 03^. Small 
deviations between the Id energies calculated using the 
renormalized Id coupling constant gid and the 3d ener- 
gies calculated using the short-range potential V SR are 
visible; we attribute these to the finite range ofV SR . The 
deviations should decrease with decreasing range tq of the 
short-range potential V SR . On the other hand, r deter- 
mines to first order the energy-dependence of the scat- 
tering length a^d- Thus, usage of an energy-dependent 
coupling constant gid [see Eq. (J5J)] should also reduce the 
deviations between the Id energies and the 3d energies 



calculated using the short-range potential V s R Q. Such 
an approach is, however, beyond the scope of this paper. 

We conclude that the renormalization of the effective 
Id coupling constant gid and of the Id scattering length 
aid are crucial to reproduce the results of the 3d Hamil- 
tonian H$d when a 3 d > a p and when a 3 d is negative. 

In addition to treating the Id many-body Hamilto- 
nian using the (FN-)DMC technique, we solve the Id 
Schrodinger equation using the LL equation of state 
(gid > 0) and the hard-rod equation of state (gid < 0) 
within the LDA (see Sec. IIII|) . These treatments are ex- 
pected to be good when the size of the cloud is much 
larger than the harmonic oscillator length a z , where 
a z = yJhfmuT z , that is, when a^d is large and positive 
or when a 3 d is negative. 

Dotted lines in Fig. show the Id energy per particle 
calculated within the LDA for g® d (using the LL equa- 
tion of state), while solid lines show the Id energy per 
particle calculated within the LDA for gi d , Eq. © (using 
the LL equation of state for gi d > 0, and the hard-rod 
equation of state for gid < 0). Remarkably, the LDA en- 
ergies nearly coincide with the Id many-body DMC en- 
ergies calculated using the unrenormalized coupling con- 
stant (pluses) and the renormalizcd coupling constant 
(squares), respectively. Finite-size effects play a minor 
role only for a^d *C a p . Our calculations thus establish 
that a simple treatment, i.e., a hard-rod equation of state 
treated within the LDA, describes inhomogeneous quasi- 
ld Bose gases with negative coupling constant gid well 
over a wide range of 3d scattering lengths a 3 d- 

For a^d —* 0~, that is, for large aid, the hard-rod equa- 
tion of state treated within the LDA, cannot properly de- 
scribe trapped quasi-ld Bose gases, which are expected to 
become unstable against formation of cluster-like many- 
body bound states for aid ~ l/ n i<2- Thus, Sec. [^inves- 
tigates the regime with negative 03^ in more detail within 
a many-body framework. 



VI. STABILITY OF 
QUASI-ONE-DIMENSIONAL BOSE GASES 

This section discusses the stability of inhomogeneous 
quasi-ld Bose gases with negative gid, that is, with 
03d > a% d and a 3 d < 0, against cluster formation. Sec- 
tion IVBl shows that the (FN-)DMC results for the Id 
Hamiltonian, Eq. I|2U|I . are in very good agreement with 
the FN-DMC results for the 3d Hamiltonian. Hence, we 
carry our analysis out within the Id model Hamiltonian, 
Eq. I|20p: we believe that our final conclusions also hold 
for the 3d Hamiltonian, Eq. H19[). For the inhomogeneous 
Id Hamiltonian Hid, Eq. JSUJ, the lowest-lying gas-like 
state is a highly-excited state (see Sec. Illlfl . We now 
address the question whether this state is stable quanti- 
tatively using the VMC method. 

We solve the Id many-body Schrodinger equation for 
the Hamiltonian H id , Eq. lf2*U|l . by the VMC method us- 
ing the trial wave function ipx given by Eqs. (|39|l and 
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FIG. 7: VMC energy per particle, Eid/N~huj p , as a function 
of the variational parameter a z for TV = 5, A = 0.01 and 
aid/o, p — 1.0326 (pluses), 2 (asterisks), 3 (diamonds) and 4 
(triangles). An energy barrier is present for aid/a p — 1.0326 
and 2, but not for a ld /a p = 4. 

(140(1 . This many-body wave function has the same nodal 
constraint as a system of TV hard-rods of size aid- How- 
ever, contrary to hard-rods, for interparticle distances 
smaller than aid the amplitude of the wave function in- 
creases as \z\ decreases. This effect arises from the at- 
tractive nature of the Id effective potential and gives rise 
in the many-body framework to the formation of cluster- 
like bound states as the average interparticle distance is 
reduced below a certain critical value. 

Figure shows the resulting VMC energy per particle, 
Eid/N— frujp, for TV = 5 and A = 0.01 as a function of the 
Gaussian width a z for four different values of aid- For 
aid/dp = 1-0326 and 2, Fig. shows a local minimum 
at ctz^min ~ a z . The minimum VMC energy nearly coin- 
cides with the FN-DMC energy (see also Fig. [SJ, which 
suggests that our variational wave function provides a 
highly accurate description of the quasi-ld many-body 
system. The energy barrier at a z ss 0.2a 2 decreases with 
increasing aid, and disappears for aid/a p w 3. We in- 
terpret this vanishing of the energy barrier as an indica- 
tion of instability [31j. For small a±d, the energy barrier 
separates the lowest-lying gas-like state from cluster-like 
bound states. Hence, the gas-like state is stable against 
cluster formation. For larger aid, this energy barrier dis- 
appears and the gas-like state becomes unstable against 
cluster formation. 

We stress that our stability analysis should not be con- 
fused with that carried out for attractive inhomogeneous 
3d systems at the level of mean-field Gross-Pitaevskii the- 
ory j33. In fact, a mean- field type analysis of inhomoge- 
neous Id Bose gases does not predict stability of gas-like 
states [3i| . In our analysis, the emergence of local energy 
minima in configuration space is due to the structure of 
the two-body correlation factor ji(z) entering the VMC 



FIG. 8: VMC energy per particle, E\dj 'TV ~T%u) p , as a function 
of the variational parameter a z for a\d/a p — 1.0326 (cor- 
responding to the unitary regime), A = 0.01, and N — 5 
(pluses), 10 (asterisks), and 20 (diamonds). (The TV = 5 data 
are also shown in Fig. 0) The height of the energy barrier 
decreases with increasing TV. Horizontal solid lines show the 
corresponding energies for TV = 5, 10 and 20 obtained us- 
ing the FN-DMC technique, which are in excellent agreement 
with the VMC energy obtained for a z — a Z:m i„. 

trial wave function tpx, Eqs. (|39|l and (|40|) . It is a many- 
body effect that cannot be described within a mean-field 
Gross-Pitaevskii framework. 

To additionally investigate the dependence of stability 
on the number of particles, Fig.[S]shows the VMC energy 
for A = 0.01 as a function of the variational parameter 
a z for different values of TV, TV = 5, 10 and 20. The scat- 
tering length aid is fixed at the value corresponding to 
the unitary regime, aid = 1.0326a p . Figure [S] shows that 
the height of the energy barrier decreases for increasing 
V. Figures and |H1 suggest that the stability of Id Bose 
gases depends on aid and V. To extract a functional 
dependence, we additionally perform variational calcu- 
lations for larger TV and different values of A and aid- 
We find that the onset of instability of the lowest-lying 
gas-like state can be described by the following criticality 
condition 

v^— ~ 0.78 , (46) 
a p 

or, equivalently, by ^/~Naid/a z ~ 0.78. Our Id many- 
body calculations thus sugg est that the lowest-lying gas- 
like state is stable if V N\aid/a p < 0.78, and that it is 
unstable if \J N\aid/a p > 0.78. The stability condition, 
Eq. (|46|) , implies that reducing the anisotropy parameter 
A should allow stabilization of relatively large quasi-ld 
Bose gases. 

To express the stability condition, Eq. I|46|l . in terms of 
the Id gas parameter nid«id, where nid denotes the linear 
density at the trap center, we approximate the density 
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FIG. 9: TG density [Eq. \Y2('tt . dashed lines] as a function of z 
together with VMC density (solid lines), obtained by solving 
the Id many-body Schrodinger equation, Eq. 1221 . for N = 5 
and aid/dp — 3.6, for N — 10 and aid/a p = 2.6, and for 
TV = 20 and a 1 d/a p = 1.8. The TG density at the center of 
the trap, z = 0, deviates from the VMC density at the center 
of the trap by less than 10 %. 

for negative gid by the TG density, Eq. Figure 

compares the TG density with that obtained from the 
VMC calculations for N = 5, 10 and 20 and values of 
aw/ o-p close to the criticality condition, Eq. (|4^|l . The 
density at the center of the trap is described by the TG 
density to within 10 %. Sin ce the TG density at the 
trap center is given by y2N/{%a z ) [see Eq. the 
stability condition, Eq. I|46[l . expressed in terms of the 
Id gas parameter reads ri\da,\d ^ 0.35. 

VII. CONCLUSIONS 

This paper presents a thorough study of the properties 
of inhomogeneous, harmonically-confined quasi-ld Bosc 
gases as a function of the 3d scattering length a^d- The 
behavior of confined Bose gases strongly depends on the 
ratio of the harmonic oscillator length in the tight trans- 
verse direction, a p , to the interaction range r$ and to the 
average interparticle distance 1/n 1 / 3 , where n denotes 
the 3d central density. 

Quasi-ld bosonic gases have been realized experimen- 
tally in highly-elongated harmonic traps. The strength of 
atom-atom interactions can be varied over a wide range 
by tuning the value of the 3d s-wave scattering length 
0,3d through application of an external magnetic field in 
the proximity of a Feshbach resonance. For vq <C a p , the 
scattering length 03d determines to a good approximation 
the effective Id scattering length aid an d the effective Id 
coupling constant g\d, which can be, just as the 3d cou- 
pling constant, tuned to essentially any value including 
zero and ±00. By exploiting Feshbach resonance tech- 



niques, one should be able to achieve strongly-correlated 
quasi-ld systems. The strong coupling regime is achieved 
for 1/n 1 / 3 ^> a p , it includes the TG gas, where a system 
of interacting bosons behaves as if it consisted of non- 
interacting spinless fermions, and the so-called unitary 
regime, where the properties of the gas become indepen- 
dent of the actual value of a^d- In the unitary regime, 
the gas is dilute, that is, nr^ <C 1, and at the same time 
strongly-correlated, that is, n|ci3d| 3 3> 1. 

The present analysis is carried out within various the- 
oretical frameworks. We obtain the 3d energetics of the 
lowest-lying gas-like state of the system using a micro- 
scopic FN-DMC approach, which accounts for all degrees 
of freedom explicitly. The resulting energetics are then 
used to benchmark our Id calculations. Full microscopic 
Id calculations for contact interactions with renormalized 
coupling constant g\d result in energies that are in excel- 
lent agreement with the full 3d energies. This agreement 
implies that a properly chosen many-body Id Hamilto- 
nian describes quasi-ld Bose gases well. 

We also consider the LL and the hard-rod equation 
of state of a Id system treated within the LDA. These 
approaches provide a good description of the energy of 
the lowest-lying gas-like state for as few as five or ten 
particles. Finite size effects arc to a good approxima- 
tion negligible. Our detailed microscopic studies suggest 
that these LDA treatments provide a good description 
of quasi-ld Bose gases. In particular, we suggest a sim- 
ple treatment of Id systems with negative g±d using the 
hard-rod equation of state. 

Finally, we address the question of whether the lowest- 
lying gas-like state of inhomogeneous quasi-ld Bose gases 
is actually stable. We find, utilizing a variational Id 
many-body framework, that the lowest-lying gas-like 
state is stable for negative coupling constants, up to a 
minimum critical value of \gid\- Our numerical results 
suggest that the stability condition can be expressed as 
niddid — 0.35. Since our conclusions are derived from 
variational Id calculations, more thorough microscopic 
calculations are needed to confirm our findings. We be- 
lieve, however, that our findings will hold even in a 3d 
framework or when three-body recombination effects are 
included explicitly. 

While our study was performed for inhomogeneous 
quasi-ld Bose gases, many findings also apply to homoge- 
neous quasi-ld Bose gases. Furthermore, the Fermi-Bosc 
mapping II 111 Hi, which allows one to map an inter- 
acting Id gas of spin-polarized fermions to an interacting 
Id gas of spin-polarized bosons, suggests that many of 
the results presented here for quasi-ld Bose gases may 
directly apply to quasi-ld Fermi gases. 
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